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Abstract 

Analytical arguments suggest that a large class of scalar field potentials permit the 
existence of oscillons - pseudo-stable, non-topological solitons - in three spatial di- 
mensions. In this paper we numerically explore oscillon solutions in three dimensions. 
We confirm the existence of these field configurations as solutions to the Klein- Gorden 
equation in an expanding background, and verify the predictions of Amin and Shi- 
rokoff [1] for the characteristics of individual oscillons for their model. Further, we 
demonstrate that significant numbers of oscillons can be generated via fragmentation 
of the inflaton condensate, consistent with the analysis of Amin [2]. These emergent 
oscillons can easily dominate the post-inflationary universe. Finally, both analytic 
and numerical results suggest that oscillons are stable on timescales longer than the 
post-inflationary Hubble time. Consequently, the post-inflationary universe can con- 
tain an effective matter-dominated phase, during which it is dominated by localized 
concentrations of scalar field matter. 
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1 Introduction 



Oscillons are long-lived, oscillatory, spatially localized solutions that exist in a large class of 
real- valued, scalar field models with nonlinear interactions [H|2l[3j|lJ[5j|6j[7J|8j[9l 
□31 El H3 ESI EH EE! HSl 113 I2H S3, and can arise from generic initial field configuations 
[231 12H 1251 121 [2H EH I2SI I2S1 EOJ- Consequently, it is natural to ask whether oscillons are 
found in the post-inflationary universe, which is typically dominated by real-valued field(s) 
with nonlinear interactions. This question has two parts. Firstly, we need to confirm whether 
oscillonic solutions actually exist for potentials corresponding to inflationary scenarios. Sec- 
ondly, given the existence of these solutions, we can ask if the post-inflationary dynamics 
naturally leads to the generation of oscillons. In this paper, we numerically confirm previous 
analytic estimates of oscillon properties, showing that both questions can be answered in the 



From the cosmological perspective, it has long been known that the post-inflationary 
universe can host a vast array of strongly nonlinear processes, thanks to parametric resonance 
and (p)reheating [3H E21 E31 EH ESI 1^1 EH EHl [Ml SOI SU S2] - Resonance fragments the 
inflaton condensate, rendering the universe strongly inhomogeneous on sub-horizon scales, 
and our work will make it clear oscillon formation is one of a long list of phenomena associated 
with resonance (also see [21 [23 ES E3 [30]). Resonance is usually linked to reheating and 
thermalization, but oscillons are effectively large, long-lived "blobs" of scalar field matter, 
and thus provide a mechanism by which resonance can be followed by an effective matter- 
dominated phase. 

Oscillons are not exactly stable, as they are not associated with a strictly conserved 
charge in the underlying theory. In practice, however, they decay on timescales much longer 
than m _1 , where m is the mass in the scalar field Lagrangian [4"3l [HJ I2Q1 [221 122]. How- 
ever, for oscillons to have significant cosmological consequences their lifetime should also 
exceed the Hubble time at formation (H" 1 ), an interval which is also typically much longer 
than m~ l . Interestingly, while oscillonic solutions have long been known, these solutions 
typically have relatively small amplitudes and suffer from a collapse instability against long- 
wavelength perturbations, with the result that their typical lifetime t osc falls in the range 
mT x < t osc < H~ l . Recently, Amin and Shirokoff [1J described a new class of oscillon so- 
lutions, corresponding to "flat-topped" configurations with large masses, and these objects 
are stable against long- wavelength perturbations. Further, in [2], Amin showed that these 
oscillon solutions can be found in potentials that also undergo post-inflationary parametric 
resonance, which drives the formation of oscillons. 

Given that the mass and lifetime of the flat-top oscillons exceeds that of their small- 
amplitude, "Gaussian" cousins [HI2], they provide a recipe for constructing oscillon-dominated 
cosmological scenarios. Our starting point is the potential 
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with A > and (X/g) 2 < 1. A necessary condition for the existence of oscillons is that 
the potential "opens out" as (p increases from zero, hence A > 0. The <£> 6 term is required 
for the flat-top solutions, and the underlying analysis makes use of an expansion in (X/g) 
[U |2]. Although the linear stability of the oscillon solutions was investigated in p], it is 
important to investigate the robustness of oscillon solutions numerically, given that this 
system is intrinsically non-linear. In addition, although [2] provided semi-analytic results for 
the post inflationary production of oscillons in three dimensions, a detailed numerical study 
of the fully non-linear evolution was performed only in 1 + 1 dimensions. Consequently, it is 
clear that we need to numerically confirm both the existence and generation of oscillons in 
three spatial dimensions before proceeding to an analysis of their consequences. 

We note that the potential considered here is not a viable model of inflation, since 
V(<p) oc y? 6 at large (p, and a sixth order slow roll potential is clearly at odds with the data. 
Rather, V(ip) is taken to be a truncated Taylor series that describes the potential at small 
field values, or (p <C m p i, and is applicable only after inflation ends. The inflationary portion 
of the potential (typically <p > m p i) is not accurately captured by the truncated expression. 
We are further assuming that the self couplings are more important than couplings to other 
fields during the era of interest. The hard lower limit on the thermalization temperature is 
set by the generation of a cosmological neutrino background and successful nucleosynthesis 
[35] EH] while inflation is typically a (near) GUT scale phenomenon. Consequently, while 
these couplings must be non-zero they can still be very small. However, it is easy to imagine 
scenarios in which couplings to other fields drain the energy from the inflaton condensate 
before oscillons form. In this paper, however, we focus on the generic, truncated potential 
and will discuss oscillon formation in realistic inflationary scenarios including the impact of 
couplings between the inflaton and other fields in a future paper. We will further restrict out 
attention to scenarios where (X/g) 2 < 10 -1 , which corresponds to insisting that the tp 6 term 
in the Taylor expansion is large. In realistic scenarios where this condition is not satisfied, 
oscillon formation is still entirely possible, but the solutions may not be well described 
by the approximately single frequency solutions found when (X/g) 2 <C 1. Finally, this 
analysis is performed in an expanding universe but ignores "local" gravitational interactions 
associated with the oscillons, which can potentially lead to primordial blackhole formation 
[3Tj or gravitational wave production [38] [39] [50] ED [521 E21 EI] , and these topics will also 
be pursued in future work. 

Oscillon-like solutions in both the post-inflationary universe and more general nonlinear 
theories have been widely studied. Oscillon formation with quasi-thermal initial conditions 
in a 1 + 1 dimensional de-Sitter universe was investigated in |25j, while oscillon-oscillon 
and oscillon-domain wall interactions in 2 + 1 dimensions are discussed in [15]. Q-balls 
[55] EH] EH EE] EH] [HOI EI] have many similarities with oscillons, but differ in that they arise 
in complex valued fields with conserved charges. Conversely, the oscillons discussed here exist 
in a real valued field with no exactly conserved charge, although we do have an adiabatic 
invariant [62]. Likewise, the formation of "axitons" in the axion field during the QCD phase 
transition was investigated in [63J and pseudo-solitons (which are not necessarily oscillons) 
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Figure 1: Sample output from a 3D oscillon simulation. Contours delineate regions where 
the local energy density exceeds the average value by a factor of 10, showing localized, 
approximately spherical, overdensities characteristic of oscillons. The comoving simulation 
volume contains a significant fraction of a single, post-inflationary Hubble volume. 

have been investigated in certain supersymmetric (SUSY) hybrid inflation models [29, 30J. 
A very detailed numerical study of Q-ball formation in gravity mediated SUSY breaking 
models is given by |64j. An analytic and numerical investigation (64 3 lattice) of inflaton 
fragmentation into oscillons was carried out in [27] [28] for a class of supersymmetric models. 
For the model described by equation 0, Amin [2] provided semi-analytic predictions for 
the post-inflationary number density of oscillons in 3 + 1 and 1 + 1 dimensions while the full 
non-linear field evolution was investigated in 1 + 1 dimensions. In these 1 + 1 dimensional 
simulations, oscillons comprised ~ 80% of the total energy density at late times. Finally, 
Gleiser and collaborators treated oscillon production in the early universe from quasi-thermal 
initial conditions, reporting a 1.2 — 3% fraction of energy density in oscillons for a double-well 
model and an SU{2) model [S3], in contrast to the fully oscillon dominated phase (~ 75%) 
we describe here. 

Numerical simulations of oscillons in 3 + 1 dimensions are challenging: oscillons are 
small, localized objects of fixed physical size. Conversely, they are generated by parametric 
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resonance at comoving scales that redshift relative to the oscillon radius as the universe 
expands, so any simulation must resolve a broad range of spatial scales, with good stability 
over a number of Hubble times. Further, the properties of oscillons depend nontrivially 
on the number of spatial dimensions, so we must work with a full 3D code, without any 
symmetry assumptions. The simulations described here were performed using two different 
codes, a modified version of Defrost [SB], and PSpectRe [S7]£] Defrost uses a finite 
differencing scheme, whereas PSpectRe is a pseudo-spectral code in which the fundamental 
degrees of freedom are the momentum modes of scalar field(s), and the two codes agree well. 
A representative output from one of our simulations showing the oscillon distribution well 
after transients associated with the initial formation phase is shown in Figure [Tj 

This paper is organized as follows. Oscillon solutions and their stability are discussed in 
Section [2] In Section [3] we review the amplification of zero point fluctuations in the inflaton 
field by parametric resonance, and estimate the resulting number density of oscillons. Section 
[4] contains the main new results in this paper. There, we provide numerical results and 
compares them to the theoretical estimates. Finally, we conclude with a summary of our 
results and possible future directions. 



2 Oscillons: Analytic Expectations 



We begin by summarizing the analytic expectations for oscillon solutions in three dimensions 
as described in [1] before numerically verifying these results. The inflaton action is 
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where Q is the determinant of the metric, R is the Ricci scalar, m p i is the reduced Planck 
mass and h — c — 1. The potential V(y?) is given by equation ([I]) near its minimum: 
m is the effective mass of the inflaton at the bottom of the potential, while A and g are 
dimensionless parameters. At large values of (p > m p i, the potential must be relatively flat 
in order to support inflation. We do not make a direct connection to detailed inflationary 
phenomenology at this point. 



The field ip obeys the usual Klein-Gordon equation, 

Dtp = V{<p). 

Assuming spherical symmetry and ignoring expansion (for now) gives 

d?f — d 2 r p> — -d r (p + m 2 ip — Ay? 3 + -^-y? 5 = . 
r m z 



(3) 



(4) 



*We added support for a sixth order single field potential to Defrost, and adapted it to run in an MPI- 
based cluster environment. PSpectRe and the modified Defrost are available for download at |http : | 
//easther .physics .yale . edu/downloads .html The visualizations of our simulations are available at 
http : //www.mit . edu/~mamin/ oscillons .html and http : //easther .physics . yale . edu/f ields .html 
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Figure 2: The left panel shows the theoretically calculated radial field profiles of oscillons 
from Amin and Shirokoff pQ. The right panel shows the radial cross-section of the energy 
density of two oscillons from a fiducial run showing the flat topped, as well the Gaussian 
oscillon. The white surface represents our one parameter fit based on the theoretically 
calculated profiles. 
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Spatially localized solutions generically require a potential that is "shallower" than quadratic 
near its minimum. Given equation (jTl), this implies A > and we assume that (X/g) 2 <C 1. 
This latter constraint is not required, but is the basis of the semi-analytic treatment of the 
flat-top solutions pp. For further convenience, we take A,j < 1 and m/m p \ <C 1. The 
fiducial values we have in mind are m/m p \ ~ 5 x 10~ 6 and X,g ~ 1CT 6 , while typical field 
values satisfy tp < m/\fX. 

We are interested in field configurations that are spatially localized and periodic in time, 
and start with the following ansatz for a one parameter family of solutions indexed by a is 



. . m 



$ a (g) cosut + 0[(X/g) 



31 



(5) 



where we define u 2 = m? — (X/g) 2 a 2 and g = (X/gm)r pp. The expansion in (X/g) is 
similar to the often used e expansion (e.g. [IB]), but can access the entire allowed range of 
harmonic, spatially localized solutions. Requiring that <f{t,r) is both periodic and bounded 
can be shown to lead to the requirement that the spatial envelope satisfies 

d 2 ^ a + -d e <S> a - a 2 $ Q + hi - h 5 a = , (6) 

which is the equation for the oscillon profile pQ. 

Assuming smoothness at the origin and spatial localization provides boundary conditions 
at r = and r — > oo. The profile depends smoothly on a, which must lie between < a < a c 
for single-frequency, localized solutions where a c = a/27/160. As shown in Figure 2j a <C a c 
gives small-amplitude Gaussian profiles, while we find wide, "flat-topped" profiles as a — > a c 
PQ. At a C} the solution becomes homogenous in space, with amplitude $ ac = a/9/10. 

The stability of oscillons is determined by the slope of the energy of the oscillon vs. a 2 
curve pQ (also see [HIES])- The presence of the ip G in the potential allows this slope to change 
sign at a s w 0.15. In the a<a s limit (small-amplitude "Gaussian" profile), 3 + 1 dimen- 
sional oscillons suffer from a collapse instability when the wavelength of the perturbations 
is comparable to the width of the oscillons |1|. This instability disappears for a s < a < a c . 
Consequently, we do not expect to see persistent solutions with a < a s . Conversely, in the 
a — > a c limit, the oscillons are extremely wide, and efficiently transfer energy to (k ~ m) 
perturbations [HI22]- Wide oscillons are locally analogous to a homogeneous, oscillating field 
and can pump high frequency resonant modes (as described in the next section). As a result, 
we expect an upper limit on the width of of the oscillons produced in our simulations. 

Extracting the profile $ Q from numerical simulations is a nontrivial task, since the oscil- 
lons can be somewhat asymmetric and the field profile is modulated by the time dependent 
cos(o;t) term. Conversely, the energy-density profiles, p a does not vary significantly with 
time, and is related to $ a by 
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The two critical values of p a , corresponding to a c and a s , are 



p c = m 4 (9A/20^ 2 ) [1 + 0[(\/g) 2 } , 

m 4 A (8) 



0.23- 



9 2 



We expect to see long lived oscillons in the regime p ag < p a < p ac . 

Given that we are assuming a radial profile, this analysis gives a match between the 
width (or, more conveniently, the cube root of the volume) and core density of the oscillons 
seen in our simulations, and this will be the relationship we seek to confirm numerically. 
For convenience, we take individual oscillons to occupy the region where the energy density 
is greater that 1/e of its value at the peak. Unlike most cases discussed in the literature, 
the width is a non-monotonic function of the core amplitude, with the increasing width at 
large core amplitudes being driven by the <p 6 term [1]. The relationship between the height 
and width of the oscillon energy densities is plotted in Figure [71a). Furthermore, beyond p c 
the width is a two-valued function of the amplitude, but still approaches the homogeneous 
solution via the flat-top profiles. 

When the expanding background is included, the solutions change character at r* ~ 
aKg/X)!!]^ 1 , becoming oscillatory in space [U [25]. Thus, there is an upper limit on the 
oscillon width in an expanding universe, even though equation ^ has solutions of arbitrary 
width. Physically, oscillons with width > r* are stretched by the expansion while the flat 
space treatment is reliable when the oscillon width obeys r e < r*. Finally, the expanding 
background slowly extracts energy from the oscillons in the form of outgoing radiation [25j 
|9j [1], in addition to the tiny, radiative loss not captured by a perturbative analysis [33]. 



3 Resonance and Oscillon Generation 

The formal existence of oscillon solutions is not sufficient to ensure that oscillons are actually 
generated in the very early universe. Immediately following inflation, the universe is homo- 
geneous, with the inflaton oscillating about the minimum of V(ip). However, the same higher 
order terms in the potential that support the oscillon configurations also drive parametric 
resonance (see for example [SU E21 El [Ml E51 ESI EH ESI EHl SOI SD S21 12])- Consequently, 
field perturbations at wavelengths within the resonance bands grow quasi-exponentially. For 
the potential here, parametric resonance in an expanding background and its consequences 
for oscillon formation were discussed by Amin in [2J, along with an approximate prediction 
for the resulting number density. We review the relevant results below before testing them 
numerically. 
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3.1 Background and Resonance and in the Linear Regime 



While resonance is often studied in multifield models, it also occurs in single field models 
with non-quadratic potentials. The resonance bands for the potential in equation are 
described in [2j. The homogeneous field obeys 



dt<p + 3Hd t ip + V'{(p) = 0, 



H 



3^pi l 



-d 2 t (p + V{0) 



(9) 



where H is related to the scale factor a in the usual way, H = a(t)/a(t). At the bottom of 
the potential, we have 



(fit) 



H 



Hi 



cos(uot) 



(10) 



a 3 (t) 



and oj 2 ~ (m 2 — ^<^ 2 + . . .) when (p is small enough that the potential is dominated by the 
quadratic term. We define a(tj) = a» = l,<p(ti) = (pi, H(tj) = Hi and choose (pi (and hence 
ti) to be the value of the homogenous field at which modes first enter the resonance band 
responsible for oscillon formation (see [2] for details), 



<Pi 



3A 



-m . 



'IV, 



During coherent oscillations, the universe is briefly potential-dominated during each cycle, 
with a instantaneously positive. Consequently, we chose the initial phase by setting a(tj) = 0, 
or 

m 




m. 
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3.2 Initial Conditions 



Given our choice of t{ (via equation (11)) we begin our simulations well after the end of 
inflation, and we must provide initial conditions for the field fluctuations, 5ip(ti). We express 
these via a Fourier amplitude 8<fk, where k is the physical wavenumber at ti. Perturbations 
that leave the horizon during inflation are substantially modified relative to their original 
Bunch-Davies form, and evolve further following horizon re-entry after inflation [68, 69J. 
However, oscillon formation is dominated by modes that never left the horizon, and we set 
the initial perturbation spectrum consistent with the zero point fluctuations of the field. 
Treating these fluctuations as a classical, Gaussian, random field, a typical fluctuation at a 
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lengthscale k 1 is given by (ignoring interaction terms): 

kni^kimY 2 - ^= (13) 

where uj\ = k 2 +m 2 + 0[H 2 ] (see [70], for example), and H is the Hubble parameter. Since we 
assume that Hi <C m, we will ignore the Hi piece. We treat the fluctuations classically at all 
times, even though their occupation number is initially zero. However, parametric resonance 
amounts to particle production, and occupation numbers grow rapidly once resonance begins, 
justifying our use of classical field theory for these simulations [ID] . 



Our potential also has narrow resonance bands for (p > (pi (t < tj) that can change 
the initial conditions described above when k >m. We evolved the linearized perturbations 
numerically for t < ti for some "sample" cases. Some modes with k >m are indeed amplified 
but this does not appear to have a significant effect on oscillon production or evolution. We 
also checked that our numerical simulations and qualitative results are robust against large 



changes (~ 10 2 ) to the initial fluctuation amplitude. From now on we use (13) to specify 
initial conditions at t,-, and concentrate of the evolution of fluctuations after 



3.3 Resonance 

In Fourier space, the linearized field fluctuations obey 

d 2 5<p k + 3Hd t S Vk + ( J^y + V"(fp)\ 5<p k = F k [Sg^] , (14) 

where the right hand side is due to fluctuations of the metric, and the solution depends 
strongly on the relationship between k,H and V"((p). We assume V"(0) — m 2 ^> F [50^], 
and thus ignore the gravitational effects. We have checked that this assumption is valid over 
the length and time scales of interest for oscillon formation and for the parameter ranges 
considered here. The effects of field and gravitational instability in a related system in the 
post recombination universe was also discussed in [7T] . 

The instability bands and Floquet indices for the oscillon potential are derived in [2], 
and plotted in Figure [3]^a). The colors represent the real part of the Floquet exponent 
as a function of the field amplitude (p and physical wavenumber k p , in flat space. In an 
expanding universe, (p ~ cpia~ 3 ^ 2 and k p = kaT 1 . As a result, a fluctuation with a given 
co-moving momentum samples different regions of this instability chart. Given that the 
expansion rate is slow compared to m, the amplification of a given comoving mode can be 
approximated by integrating the Floquet exponent along the thin black curves shown in 
Figure |3](a). The band of wavenumbers that leads to oscillon formation is 

lli < k < 0.65 ( - j m. (15) 
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Figure 3: In Figure (a), the resonance band for the potential V((p) = m 2 (p 2 /2 — A<^ 4 /4 + 
g 2 tp 6 /6m 2 is shown. The shaded region has a positive Floquet exponent. The thin, black 
lines show the "path" following by modes as the universe expands. Our simulations begin 
at the thin, dashed line. The thick, dashed indicates the homogenous field value where the 
evolution becomes non-linear, with fluctuations becoming comparable to the homogeneous 
field. In Figure (b), the fluctuation power spectrum is shown at the point where linear 
analysis starts breaking down. The horizontal axis in Figure (a) is physical momentum, 
whereas in Figure (b) it is the co-moving momentum normalized at thin, dashed line. The 
momenta and field values are expressed in units of the mass, m. Figure (a) is adapted from 

m 
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The amplification undergone by these modes can be approximated by 



where 




Significant amplification requires /3 ^> 1, which provides an important constraint on the 
parameters X,g and m [2]. The parameter ~ fi/H where [i is the Floquet exponent 
characterizing the growth rate of fluctuations in flat space and H is the Hubble parame- 
ter. Hence the condition j3 ^> 1 follows from noting that to have rapid amplification of 
fluctuations, their growth rate /i must be larger that the Hubble rate H. Typical resonant 
modes have k < m(X/g), and ji ~ 0.1m(X/g) 2 . We have assumed (X/g) 2 ~ 0.1, while 
H ~ m(X/g) 2 /f3 ~ lCT 3 m. For these values /3 ~ 10 2 . The numerically-evaluated spectrum 
of fluctuations (linear regime) is shown in Figure |3^b). The spectrum is in excellent agree- 
ment with the Floquet analysis, justifying the approximations made in deriving the above 
expressions [2]. 

We have concentrated on the first instability band and largely ignored the higher order 
bands (for example at k p ~ \/?>m). The higher order bands are very narrow (AA; P < (X/g) 2 ) 
so modes redshift rapidly through these bands, with their maximum Floquet exponent scaling 
as /i oc (X/g) 2 . Moreover, oscillons are significantly larger than m _1 so we do not expect 
these high-fcp modes to significantly affect oscillon formation. Nevertheless, if modes pass 
through a large number of such bands or spend a long time in one of these narrow bands 
(that is, H is very small) large amplification is possible and could lead to oscillon formation. 
We checked numerically that these high-fc modes have little impact on oscillon formation for 
this scenario Q 



3.4 Number density of oscillons 

The linear analysis breaks down when dip ~ fc 3//2 (|5</9fc| 2 ) 1//2 is comparable to (p as we no longer 
have a homogeneous pump field to parametrically amplify the fluctuations. In addition, as 
Sip approaches V"{(p)/V"'{ip), we have to take into account the interaction between different 
k modes. The modes start interacting strongly with each other, ultimately forming oscillons. 
This process is difficult to follow analytically, but we can estimate the initial number density 



^We show spectra in Figure 11 - the small bump at k p ~ m 1 (for a < 2) is produced by the the 
k p ~ \/?>m band. We discuss instability further at the end of Section [4] 
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of oscillons, and its dependence on the free parameters in the potential ([!]) based on the 
linear analysis of growth of fluctuations [2]. 

As parametric resonance ceases to be efficient, there is a characteristic comoving length 
scale, fc n i where the fluctuations in the field are most nonlinear, namely the mode where the 
power spectrum peaks (see Figure [3| 
oscillons can be estimated as 



and Figure 11). Following j2], the number density of 



-a 3 ~ 



where the characteristic scale k n \ is given by 



2tt 
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with f3 = y/X(X/g)(m pl /m). Thus 



n osc ,a ~ 



-1/5 




{\/g)m 



- 3 / 5 M x3 

A m 
~g2^ 



(19) 



(20) 



Note that this is a rough estimate. Heuristically, we are merely counting the number of large 
peaks in the energy density. The linear analysis tells us how many such peaks we should 
expect in a given volume. 



4 Numerical Simulations in 3 Dimensions 

Numerically simulating the formation and evolution of oscillons in three dimensions requires 
the same tools and techniques developed during the analysis of parametric resonance and 
preheating, which greatly facilitates our work. That said, oscillon dynamics present signif- 
icant numerical challenges as the oscillons are relatively small and well separated in space, 
so we need to simulate a large volume with good spatial resolution. Moreover, we want to 
ensure that our results are cosmologically stable, and thus follow the evolution for several 
Hubble times (H~ l ~ 10 3 m _1 ), which is a much larger interval than the natural oscillation 
period of the fields (T ~ m _1 ). 



4.1 Field Evolution Programs 

To ensure that our results are independent of our numerical algorithms, we performed sim- 
ulations with a variant of Defrost [66], a finite differencing code, and PSpectRe, which 

^Note that compared to the formula presented in [2J, there is an extra factor of \/2/3 multiplying /3 
because we are starting with a different Hubble parameter Hi. 
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Figure 4: The orange curves represent the evolution of the energy fraction in oscillons as 
we vary the initial boxsize L = (200,400,600m -1 ) with iV = 256. The blue curve shows 
the same quantity when N = 384 with L = 400m -1 . Finally the black points are based 
on Defrost output at L — 400m -1 and iV = 1024. The agreement between the curves 
provides evidence for the robustness of our results. Note that this energy fraction is an 
underestimate since here we have defined the oscillon energy as the energy within a radius 
where the amplitude falls by 1/e of its core value. 



13 



uses a pseudo-spectral scheme \§T\ ■ For the cases we compared, both codes yield similar 
results and, unless otherwise specified, all results presented here were obtained with PSpec- 
tRe. We extended Defrost to include the sixth order scalar field potential and adapted 
it for use on an MPI-based cluster, running at simulation volumes of up to 1024 3 points. 
PSpectRe runs on multiple cores on a single shared-memory node, and these simulations 
were typically performed at 256 3 . 



The relative performance of both codes as a function of grid-size is discussed in 
PSpectRe typically reaches the "continuum" limit for smaller values of N than a finite dif- 
ferencing code, offsetting the numerical overhead induced by the repeated Fourier transforms 
needed to construct the nonlinear terms in momentum space. For both codes, simulations 
could require several days of wallclock time. At 1024 3 , Defrost parallelized efficiently over 
~ 12 nodes, or ~ 100 MPI threads. For some parameter values the thin, higher order reso- 
nance bands generate short-wavelength inhomogeneity after oscillon formation (see Section 



4.5). In these cases, the relevant modes must be several times larger than the Defrost grid 



resolution in order to ensure numerical fidelity. 

We ran numerical simulations using both PSpectRe and Defrost for a wide variety 
of free parameters, as well as for a range of choices of comoving boxsize L and numerical 
resolution, measured by iV 3 where N 3 is the number of points used by Defrost, twice 
the number of momentum modes resolved by PSpectRe. Figure [4] shows how the energy 
fraction in oscillons, time of emergence etc. is independent of intial box size L and the spatial 
resolution L/N as well as the code used. 



4.2 Oscillon Detection Algorithm 

Operationally, we define an oscillon to be any contiguous region of space whose energy 
density is at least 4 times the mean density. We detect oscillons by labeling each point over 
the threshold (on a position-space grid) with a unique, positive integer. We then iteratively 
set any pair of adjacent numbers to the lower value, until no further changes are possible. 
This creates a "mask" , and a given oscillon consists of the set of grid points with the same 
numerical label. Once the oscillons are detected they can be trivially counted. 

For comparison with analytic results for the oscillon profiles, we define the volume of the 
oscillon to be the region for which the energy density is within 1/e of the peak density, and 
then compute the effective radius R = (3V/47T) 1 / 3 . The total mass of the oscillon is given by 
the integral over the energy density within this volume. For steep-sided, flat-topped objects 
this definition captures a large fraction of the total energy density in the object, but can 
exclude the shoulders of a Gaussian oscillon. Consequently, when we calculate the overall 
fraction of the energy in oscillons we will also consider the region within l/3e of the peak 
density. 

In practice, we sometimes find non-oscillonic, isolated spikes that exceed our threshold 
energy density. Rather than raising the threshold, we can filter these objects by applying 
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a cut on the energy of the spikes, as they are typically factor of 100 or more less massive 
than genuine oscillon configurations. Conversely, while oscillons are typically well-separated 
at late times, interactions do occur, and this algorithm counts two interacting oscillons as 
a single object. Finally, as a matter of convenience, we ran the oscillon detection code on 
output files from PSpectRe or Defrost, rather than counting oscillons within these codes 
as they ran. 

4.3 Qualitative Results 

Unless otherwise stated, we work with "fiducial" parameters A = 2.8125 x 10 -6 , (A/g) 2 = 0.1 
and m = 5 x 10 _6 m p i. Figure [5] depicts the spatial density profile seen in a typical simulation, 
at six different times. Surfaces are drawn at a overdensity of 1, or an absolute density of 
2, in units where p = 1. Initially, over-densities form on large comoving scales <C m^ 1 , 
and collapse quickly into roughly tubular regions, which then fragment into the localized 
overdensities that we will identify as oscillons. These overdense regions are well separated, 
and slowly moving with respect to the comoving background. We see occasional interactions 
between these regions, but these are infrequent given their small spatial size (relative to the 
overall volume) and low velocity. Further, as the space expands, overdense regions occupy a 
progressively smaller comoving volume, consistent with analytic results suggesting that the 
oscillons have approximately-fixed physical size. 

Figure [6] (left panel) shows that the expansion rate of the universe is entirely consis- 
tent with a pressureless fluid, for which a(t) oc t 2//3 , during the initial phase of coherent 
oscillations, parametric resonance, and after the onset of nonlinearity. This is contrast to 
many resonant scenarios [HH E21 E31 EH [351 ESI OZl [38l [39l SQl SH S2] in which the universe 
is radiation-dominated or in some intermediate state [66] following resonance. This also 
confirms that while the potential contains terms beyond quadratic order to support the ex- 
istence of oscillons, it is dominated by the m 2 (p 2 term during the coherent oscillations phase. 
Following resonance, Figure [6] (right panel) shows that the spatially-averaged density p falls 
rapidly, while the maximum density stays roughly constant. The maximum density is found 
at the "core" of the oscillons, and is expected to be constant as the universe expands. Theo- 
retically, we do expect oscillons to decay slowly, but at rates which are negligible on the time 
scale of our simulations. To summarize, we have confirmed that the resonant phase matches 
our expectations from the linearized theory [2] , as discussed in Section [3j while the "global" 
properties of the nonlinear phase are consistent with an oscillon dominated universe, with 
localized overdensities which are a factor of 100 times larger than the average density. 

We now seek to confirm that the overdensities we see are, in fact, oscillons. We show 
the relationship between energy density and width (volume) for regions flagged as oscillons 
in the left panel of Figure [7| and see an excellent match to the theoretical expectations of 
Amin and Shirokoff [1J. Note that at late times, almost all putative oscillons have p > p s , 
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a=1.85 a=1.93 




Figure 5: Density contours (p = p) for a 1/8 sample of a simulation with fiducial parameters 
at L = 200 and iV = 256. Frames are labelled by corresponding value of the scale factor. The 
large over-dense regions are oscillons; smaller regions are simply areas in which the density 
exceeds the average value, without being oscillonic. The frames show a fixed comoving 
volume. 
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Figure 6: The left panel shows a 3//2 as a function of time t for a "fiducial" simulation. This is 
a linear function, consistent with the post-inflationary expansion mimicking that of a matter 
dominated universe. The right panel shows the maximum energy density in the simulation 
volume p m ax, and the spatially-averaged energy density p as a function of a. The former 
is roughly constant, suggesting that the maximum density of the oscillons is fixed. The 
density contrast associated with the oscillons grows as the average density is decreased by 
expansion. The initial phase of coherent oscillations during which p max ~ p is decreasing is 
clearly visible in the second plot. The densities have been normalized to unity at a = 1. 

consistent with the stability analysis in [ljjjs] An equatorial cross-section of a flat-topped 
oscillon and a Gaussian oscillon extracted from a simulation is shown in the right panel of 
Figure [2j The white layer represents the single parameter analytic fit based on equation Q . 
Clearly, a Gaussian profile would not provide a good fit to the flat-topped oscillons. 



4.4 Statistics, oscillon number and energy-density evolution 

Figure [7] (right panel) shows a histogram of oscillon energies at late times, when aj ~ 5. 
Qualitatively, we see large number of small oscillons and a weaker peak in the distribution 
at higher energies. This interpretation is somewhat subjective, but resembles the results 
from the one dimensional case, which had two clear "generations" of oscillons [2J. Oscillons 
in the right peak (more massive) form first, and the number density of first generation 



oscillons matches well with the analytic estimate of equation (20) [2j. However, we note that 



this distinction between "generations" cannot always be made cleanly, especially for large f3 



§This is in sharp contrast with the 1 + 1 dimensional case where small amplitude oscillons are stable 
against long wavelength fluctuations. As a result the left half of the corresponding "U" shaped curve is 
equally populated 1. 2 . 
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Figure 7: The left panel shows the non-monotonic relationship between the width and core 
amplitude of the oscillon energy densities, relative to the analytic prediction in [1] (solid 
black line). The width is taken as twice the radius where the energy density falls to 1/e 
of its core value. Orange points are overdensities flagged as oscillons in our simulations at 
a ~ 5. At large core amplitude we approach a flat-topped profile. The lack of oscillons at 
(g 2 /X)p C ore < {g 2 /^)Pa s ~ 0.23 is consistent with the stability analysis of Amin and Shirokoff 
PP. The right panel shows a histogram of oscillon energies with the width r e . The parameters 
used were A ~ 2.8 x 1CT 6 , (A/g) 2 = 10 _1 and m/m p \ = 5 x 10~ 6 . All variables are expressed 
in units of appropriate powers of mass, m. 
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123456 123456 
scale factor: a scale factor: a 

Figure 8: The left panel shows the number density of oscillons as a function of scale factor. 
Oscillons emerge at a n \ ~ 1.5. Smaller oscillons are transient, leaving a roughly constant 
(comoving) number density of large-amplitude oscillons. The dashed line is the analytic 
stability estimate from Amin [2]. The energy density in oscillons is shown on the right. The 
fraction of the total energy density in oscillons, defined here as regions in which p > p CO re/3e 
for individual oscillons, is shown in the right panel. The parameters used were A ~ 2.8 x 1CT 6 
(or (3 = 10 2 ), (X/g) 2 = 0.1 and m/m p \ = 5 x 10~ 6 . Number density is expressed in units of 
m . 



values. For this paper, we do not distinguish between the two generations. 



Figure [8] (right) shows the comoving number density of oscillons as a function of time. 
Oscillons are produced at a n \ ~ 1.7, and the number density peaks there. As noted above, 
small oscillons are unstable in three dimensions, but the number density stabilizes within 
a Hubble time. The remaining large oscillons have a fixed physical size and thus shrink, 
relative to the comoving simulation volume, and their typical separation scales with the 



expansion. The dashed line is the analytic estimate based on equation (20) [2J. The fraction 
of energy density in oscillons is shown in Figure |8](left). Note that the fractional energy 
computed here underestimates the total energy associated with the oscillons. At late times 
in our simulations, the peak density can be at least 100 times larger than the average. Even 
if we include points whose energy is within l/3e of the maximum, we still ignore points 
for for which density contrast is of order 10. Qualitatively, the vast majority of the energy 
density is accounted for by the oscillons. 

We further check our expectations by varying A and g, and verifying that the numerical 
results reproduce the scalings seen in semianalytic estimates of the co-moving number density 
of large oscillons [2j. We set L = 400m" 1 and N = 256 and express the number density as 
a function of j3 = y/X(X/g)(m p i/m), since j3 roughly characterizes the strength of resonance 
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Figure 9: The above Figure shows the variation of the co-moving number density with /3 = 
\ 3 / 2 g(m p i/m) which characterizes the ratio between the maximum growth rate of fluctuations 
in flat space and the Hubble parameter. The dashed curve is the analytic estimate provided 
by Amin [2]. The orange points are for (X/g) 2 = 0.1 and the black ones are for (X/g) 2 = 0.2. 
The dependence on (X/g) is accurately captured by the analytic estimate, which results in 
the orange and black points lying on the same curve. However the dependence on the relative 
growth rate of fluctuations, (3 is stronger than in the analytic estimate. The deviation from 
the analytic estimates can be attributed to oscillon-oscillon collisions at small Hubble values 
(large f3). 



in the linear regime. Figure [9[b) shows the final number density of oscillons at late times, 
after the oscillon number density has stabilized. We show the results as a function of (3 for 
(X/g) 2 = 0.1 (orange) and (X/g) 2 = 0.2 (black), while A runs from 1.25 x 10~ 6 to 1.125 x 10" 5 . 
The dashed curve is the analytic estimate derived by Amin in [2]. Remarkably, despite the 
non-linear dynamics of oscillon formation, the analytic estimate is good to within a factor of 
a a few and correctly captures the parameter dependence. In particular, the variation with 
(X/g) is almost exactly described by the scaling of the vertical axis. 

The analytic result over-estimates the number of oscillons at large (3, where the expansion 
rate is slow. Consequently, one can infer that oscillon-oscillon interactions are more common 
at formation, which appears to lead to mergers and fewer, larger oscillons and this process 
was not part of our estimate. The outcome of a collision depends on the phase, amplitude, 
velocity as well as shape of the oscillons (see [15J). The collisions can be almost completely 
elastic ones or highly inelastic ones, though we are almost always left with one or more 
oscillons. Details of oscillon-oscillons collisions and fragmentation will be considered in 
future work. 

As noted previously, the universe is effectively oscillon-dominated after nonlinearity sets 
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Figure 10: The left panel shows the percentage of energy density in oscillons at the 
end of our simulations for different values of the relative growth rate of fluctuations: 
(5 = \f~X{X/ 'g)(m p i/m). The black and orange points correspond to simulations with param- 
eters {X/g) 2 = 0.2 and 0.1 respectively. The fraction of energy density in oscillons decreases 
at low f3 (large H) because resonance is not efficient enough to significantly amplify initial 
fluctuations. The right panel shows the mean energy per oscillon. Note that it increases 
with increasing /3 (decreasing H). This is consistent with the idea that at large /3 oscillons 
form from longer wavelength fluctuations, thus yielding fewer, larger oscillons. The orange 
and black points are aligned as expected, because of the scaling with (X/g) on the vertical 
axes. Energy is expressed in units of the mass, ra- 



in. In Figure 10 (left), we quantify this statement. At very small /3, we cannot generate 
oscillons efficiently, so the oscillon contribution initially increases with 0, but eventually 
saturates. In Figure 10 (right), we show the average energy per oscillon as a function of 
P. As (5 increases (H decreases), we get fewer, more massive oscillons. This is consistent 
with the shape of the resonance band seen in Figure [3|a): for small H the longer wavelength 
modes have a chance to get amplified before the "sweet spot" of the Floquet band is traversed 
by k ~ (X/g)m modes. Further, as discussed above, small H also allows for a larger number 
of interactions between oscillons, possibly leading to mergers and hence bigger oscillons. 



4.5 Spectrum of linear and non-linear fluctuations 

As we have seen the initial, linear growth of fluctuations is accurately and conveniently 
described in Fourier space, whereas oscillons find a natural description in position space. 
Nevertheless, it is instructive to look at the spectrum of fluctuations through the entire 



evolution. Figure 11 shows the spectrum of fluctuations for the parameters A ~ 2.8 x 10 
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Figure 11: The figures show the evolution of the spectrum of field fluctuations for A = 
2.8125 x 1CT 6 and (X/g) 2 = 0.1 (left panel) and (X/g) 2 = 0.2 (right panel). The orange and 
red curves show the linear evolution of fluctuations through resonance. This broad band is 
responsible for oscillon production. The small bump seen at physical momentum k p ~ m is 
due to higher-order resonance bands excited by the homogeneous field initially as well as by 
very wide oscillons at late times. The power increase in power at k p ~ y/3m is clearly visible 
in the right panel. 
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and (A/g) 2 = 0.1 (left panel) and (A/g) 2 = 0.2 (right panel). The orange and red curves 
represent the spectrum in the linear regime, whereas the yellow and green ones are deep in 
the non-linear regime. The biggest feature in the spectrum is from the band that produces 
oscillons initially and at late time shows the Fourier content of spatially localized oscillons. 

There is an additional, tiny bump in the spectrum at the physical momentum k p ~ m. 
A narrow, higher-order resonance band at k p ~ y/3m generates this bump during the initial 
phase of coherent oscillations. That power is redshifted and this small bump moves to the 
left. At late times, long after resonance due to the homogenous field has ended, very wide 
oscillons regenerate power at k p ~ \/3m, as they effectively behave as a homogeneous field. 
This can be seen in the slight increase in power below k p ~ \/3m (left panel). This Floquet 
instability at k p ~ Vom [TJ 122] can potentially destabilize the oscillons and drain their energy. 
However, these fluctuations are catastrophic for the oscillon only if fluctuations cannot escape 
the oscillon before being significantly amplified. This yields a relation between the rate of 
growth of fluctuations and the width of the oscillons. We refer the reader to [22] for further 
details. As seen in Figure TTJ we see this effect in k p ~ \^3m fluctuations clearly for case 
(A/g) 2 = 0.2 (right panel) but not for the case (A/g) 2 = 0.1 (left panel) at late times. This 
is expected since the Floquet exponent for the k ~ \f?>m instability band scales as oc (A/g) 2 . 
We note that for many runs we cannot resolve the k p ~ instability towards the end of 
the simulation because of the degraded resolution from the expanding background. For some 
sample higher resolution runs, we did not see a significant change in the number density of 
oscillons. Nevertheless, the effect of this instability needs to be investigated further. 



5 Discussion 



In this paper we have numerically investigated the post-inflationary emergence of oscillons 
in a class of single-field potentials in 3 + 1 dimensions. Beginning our simulations with a 
homogeneous inflaton condensate and the zero point fluctuations of the inflaton, resonance- 
driven inflaton fragmentation leads to the copious production of oscillons. In these scenarios, 
inflation is followed by a period of oscillon domination, during which the dominant component 
of the universe is localized "blobs" of scalar field matter. The universe is pressureless, leading 
to an effective matter-dominated phase following parametric resonance. The number density 
and individual characteristics of oscillons closely match the theoretical predictions of [JJ and 

This paper forms part of a wider investigation of the existence and consequences of 
oscillon-like field configurations in the very early universe. This work is based on a potential 
which is symmetric as ip — > —if, truncated at the third nontrivial term in the Taylor ex- 
pansion which is taken to be large, in order to make contact with the analytic discussion of 
[TJ Ej . Our next step in this program will be to relax these conditions, and demonstrate the 
existence of persistent, oscillon-like configurations for realistic inflationary potentials, and 
to explore the cosmological consequences of a post-inflationary, oscillon-dominated phase. 
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Looking further ahead, we wish to understand the impact of the higher-order resonance 
bands on oscillon dynamics, account for the effects of couplings to additional fields, consider 
the impact of gravitational effects on the formation and evolution of isolated oscillons, and 
the interactions between oscillons. 
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